Dark Halo Shapes and the Fate of Stellar Bars 
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Q ■ ABSTRACT 

(f) ■ We investigate the stability properties of trajectories in barred galaxies with mildly triaxial 

halos by means of Liapunov exponents. This method is perfectly suitable for time-dependent 3-D 
potentials where surfaces of sections and other simple diagnostics are not applicable. We find 
that when halos are centrally-concentrated most trajectories starting near the plane containing 
^ ■ the bar become chaotic. The spatial density distribution of these orbits does not match that of 

the bar, being overextended in- and out-of-the plane compared to the latter. Moreover, the shape 
of many of the remaining regular trajectories do not match the the bar density distribution, being 
too round. Therefore, time-independent self-consistent solutions are highly unlikely to be found. 
When the non-rotating non-axisymmetric perturbation in the potential reaches 10%, almost all 
trajectories integrated are chaotic and have large Liapunov exponents. No regular trajectories 
aligned with the bar have been found. Hence, if the evolution of the density figure is directly 
related to the characteristic timescale of orbital instability, bar dissolution would take place 
on a timescale of few dynamical times. The slowly rotating non-axisymmetric contribution to 
the potential required for the onset of widespread chaotic behavior is remarkably small. Even a 
potential axis ratio of 0.99 results in large connected chaotic regions dominating the space of initial 
conditions. Systems consisting of centrally-concentrated axisymmetric halos and stellar bars thus 
appear to be structurally unstable, and small (~ 1%) deviations from perfect axisymmetry should 
result in a bar dissolution on a timescale significantly smaller than the Hubble time. Since halos 
found in cold dark matter simulations of large scale structure are both centrally-concentrated 
and triaxial it is unlikely that stellar bars embedded in such halos would form and survive unless 
the halos are modified during the formation of the baryonic component. 

Subject headings: instabilities — stellar dynamics — galaxies: evolution — galaxies: halos — galaxies: 
structure — cosmology: dark matter 

1. Introduction and motivation 
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The main reason for this is that the breaking of 
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axial symmetry introduces gravitational torques, 
whose action can be described in terms of a non- 
local viscosity (e.g., Larson 1984; Shlosman 1991). 
This causes accelerated redistribution of mass and 
angular momentum. Disk halo interaction is also 
increased dramatically in the presence of a bar. 
This is particularly so if the halo is also triaxial 
with significant contribution to the density in the 
inner regions — that is if it is centrally concen- 
trated. Such a strong interaction raises questions 
about stability of the least massive object in this 
configuration, the stellar bar. It also has clear 
implications on disk formation and morphological 
evolution from high redshifts down to the local 
universe. In this paper, we analyze the orbital 
stability of trajectories evolving in the potential 
of a stellar bar embedded in massive halos of var- 
ious central concentrations and asymmetries and 
discuss a number of corollaries, in order to infer 
the overall stability of the bar. 

Triaxial systems must be built by dynamical 
trajectories that conserve, at least approximately, 
invariants of motion, if they are to remain in quasi- 
steady state. This is necessary if their configura- 
tion (i.e., physical) space density is to match the 
triaxial shape of the system. For example, a non- 
rotating steady state system, whose distribution 
function depends only on the energy, is necessarily 
spherical. While this is the only "global" integral 
of motion that always exists in time-independent 
systems, it is not sufficient to maintain a triax- 
ial shape. Moreover, additional global integrals of 
motion are often associated with spatial symme- 
tries, which are lacking in triaxial systems, render- 
ing the problem of finding self-consistent solutions 
highly non-trivial (for a discussion of these issues 
for the case of triaxial elliptical galaxies see, e.g., 
Merritt & Fridman 1996 and Merritt & Valluri 
1996). In some cases, symmetries different from 
simple spatial ones can exist, leading to global in- 
tegrals of motions that can maintain the triaxial 
structure. Thus is the situation with Stackel po- 
tentials (e.g., de Zeeuw 1985). In the core region 
of such potentials the effective symmetry is the 
near-homogeneity of the density distribution. The 
potential can then be approximated as a quadratic 
form where the motion is separable in Cartesian 
coordinates (in fact this type of system is not 
only separable but even linear) . Generic centrally- 
concentrated potentials do not have such symme- 



tries, the oscillations in the different degrees of 
freedom are, in general, coupled and no global in- 
tegrals of motion exist. This makes plausible a 
situation whereas most trajectories in the central 
regions conserve only energy, in which case self- 
consistent solutions become impossible. 

Centrally-concentrated mass distributions have 
been invoked by several authors in proposing a 
mechanism for secular evolution in galaxies, both 
in the case of slowly rotating elliptical galaxies 
(e.g., Norman, May & van Albada 1985; Merritt 
& Fridman 1996; Holley-Bokelmann ct al. 2002) 
and, in the context of rapidly rotating bars, in disk 
galaxies (e.g., Pfcnniger & Norman 1990; Norman, 
Sellwood & Hasan 1996). For example, it has been 
argued that barred galaxies that are not initially 
centrally-concentrated may acquire a central mass 
concentration by accreting gas into the central re- 
gion. This creates the necessary coupling between 
the degrees of freedom so as to destroy the inte- 
grals of motion for a large enough fraction of orbits 
and dissolve the bar — which is then replaced by a 
bulge-like structure. It appears, however, that the 
large central masses required for such a scenario to 
work, specifically the central black holes, are not 
confirmed by observations (in preparation). An 
analogous situation should however transpire if the 
coupling between the degrees of freedom is medi- 
ated by the existence of a centrally-concentrated 
halo dominating the inner density distribution. 
This would be the case for systems with halos of 
the type found in cosmological simulations of the 
cold dark matter (CDM) scenario of structure for- 
mation (e.g., Navarro, Frenk & White 1997; Moore 
et al. 1999). 

Halos identified in cosmological simulations are 
also invariably found to be triaxial (e.g., Warren 
ct al 1992; Cole & Lacey 1996). In general, the 
halo being built of non-dissipative material, will 
be much more slowly rotating than an embedded 
baryonic bar formed in its central region. The 
introduction of a bar, i.e., an additional triax- 
ial configuration, thus ensures that even energy 
is not conserved in any uniformly rotating frame 
of reference. The loss of an additional global (time 
translation) symmetry makes it even more likely 
that a majority of trajectories will be chaotic with 
nearly random motion, and hence will not sup- 
port the existence of the bar. Indeed earlier ex- 
ploratory work on this issue suggested that this 



2 



is the case even for non-axisymmctric halos with 
moderate constant density cores (El-Zant & Has- 
slcr 1998). This was confirmed by self-consistent 
simulations performed by Ideta & Hozumi (2000) 
for systems with centrally-concentrated axisym- 
metric halos, but with dominant disk contribution 
to the mass distribution of the inner regions. 

It is the goal of this paper to investigate the 
prospects of bar survival for different halo shapes 
in detail. We will do this by analyzing the response 
of the orbital structure supported by the system to 
various degrees of halo central concentration and 
asymmetry. Since it has been shown (Dubinski 
1994) that the settling of a baryonic component in 
a triaxial halo can significantly reduce the initial 
non-axisymmetry born of dissipationless collapse, 
we will only consider rather small departures from 
axisymmetry; deviations from unity in the poten- 
tial axis ratio of 10% or less, which would induce 
cllipticitics in galactic disks of still smaller magni- 
tude and are, therefore, consistent with observa- 
tions of even present day galaxies (for a summary 
of observational results concerning halo shapes see, 
e.g., El-Zant & Hassler 1998; Tremaine & Ostriker 
1999). 

The stability of motion along trajectories in our 
models can be quantified by a variety of methods 
developed in the dynamics of nonlinear systems 
(e.g., Hilborn 1994). In particular, we shall em- 
ploy the method of Liapunov exponents. This has 
the advantage of providing a timescale that is as- 
sociated with an instability, i.e., an e-folding time, 
is much simpler to implement with high accuracy 
than other methods (e.g., Laskar's frequency anal- 
ysis: Papaphilipou & Laskar 1998) and is suited 
for far from integrable systems where most of the 
phase space is occupied by highly chaotic trajec- 
tories — as will turn out to be the case for some 
of our models. It is readily generalized to higher 
dimensional, intrinsically time-dependent systems 
- cases where other simple diagnostics, such as 
surface of sections are not applicable. It will be 
one of our goals to see how closely the Liapunov 
timescales correlate with their time-averaged den- 
sity in configuration (i.e., physical) space. Intu- 
itively, it is expected that those trajectories that 
have small e-folding times will fill a large region of 
configuration space, generally not coincident with 
the bar density distribution. Therefore, they are 
unlikely to serve as the main contributors in build- 



ing a bar. If most, or almost all, trajectories are 
of this type then one may plausibly conclude that 
a bar is not sustainable. Sufficient conditions for 
a trajectory to belong to this group include the 
following. The volume a trajectory occupies in 
configuration space exceeds substantially that of 
the bar at the radius where the trajectory starts. 
The trajectory spends most of the time outside 
the 3-D figure of the bar, for example at vertical 
excursions larger than the bar's semi-minor axis. 
The trajectory is too isotropic in the plane, i.e. is 
too round to match the bar. 

The methods which we employ in evaluating the 
Liapunov exponents and the interpretation of their 
values are discussed in the next section, the techni- 
cal details are deferred to the appendix. The spa- 
tial density distribution of the integrated trajecto- 
ries will be examined with the aid of the configura- 
tion space grid described in Section 3. The system 
parameters of the galactic models and the units 
used here can be found in Section 4, and the initial 
conditions and integrator employed are described 
in Section 5. Section 6 contains the presentation 
of the main results of this paper. It addresses 
the stability properties of trajectories as well as 
their spatial distributions and the correspondence 
of these two properties for models with various 
halo central concentration and triaxiality. We will 
also be examining the vertical stability of trajec- 
tories and examine its dependence, as well as that 
of the general spatial structure, on the integration 
timescale. For all models we obtain the maximal 
Liapunov exponents, which is faster to calculate 
than the whole set and is usually indicative of the 
degree of instability of a given trajectory. For se- 
lected models, we calculate the full set of Liapunov 
exponents. They are indicative of conservation of 
invariants along trajectories. In the general time- 
dependent potentials studied here, all the expo- 
nents can be non-zero. Finally, in Section 7, we 
examine the correlation between the energy diffu- 
sion and the instability properties of the trajecto- 
ries as determined by the values of the Liapunov 
exponents, for one of the models. Our conclusions 
are discussed and summarized in Section 8. 
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2. Liapunov exponents and stability of 
trajectories 

Motivated by the discussion of the previous sec- 
tion, we would like to estimate the fraction of ini- 
tial conditions leading to unstable (chaotic) trajec- 
tories and the associated timescale. For this pur- 
pose we will calculate time-dependent Liapunov 
exponents for a large set of initial conditions and 
display them on grayshade diagrams. In the ap- 
pendix we describe in a more formal manner how 
the exponents are obtained, here we restrict our- 
selves to some general comments, helpful in the 
interpretation of the results of the next section. 

Liapunov exponents compare the asymptotic 
rate at which the distance in the phase space, 
|| (5X(i,X(0)) ||, between initially adjacent trajec- 
tories starting with initial conditions around X(0) 
to the exponential (see for example Lichtenberg & 
Lieberman 1995). They may therefore be written 

1. II <JX(t,X(0)) || 



lim - log ■ 

t^oo t 



II sx(o) || • (1) 

To each initial condition X(0) and perturbation 
6, there corresponds a Liapunov exponent, and to 
this a characteristic exponential timescale (its in- 
verse). Thus a mapping is created between the 
space of initial conditions and the stability associ- 
ated with them. 

If || (5X(£) ||~ t, as in the case, for example, of 
regular motion represented by action angle vari- 
ables (e.g., Binney & Tremaine 1987) — such that 
J = Ci and = tot + C2, with Ci and C2 con- 
stant vectors — the exponents converge to zero 
as - ^p. If, on the other hand, || SX(t) ||~ e\ 
as in the case of unstable systems, then the max- 
imal exponent tends to a finite limit (this limit 
can be shown to exist: see, for example, Eck- 
mann & Ruelle 1985). One important difference 
between these two cases is that exponential in- 
stability signals deviation normal to the trajecto- 
ries and, therefore, implies variation in the action 
variables and not only the phase (angle) variables. 
This can lead to time dependence. For example, 
let us suppose one of the action variables is the 
oscillations of the trajectory with respect to the 
z-axis. If this is conserved, then trajectories start- 
ing near the z = plane will remain there. If it 
is not conserved, then they will eventually venture 
out of that plane. This may have important physi- 



cal physical consequences, a system originally con- 
fined to a plane may become three dimensional. 
How fast will such evolution in the statistical prop- 
erties of trajectories take place will, in general (but 
not in all cases), depend directly on the timescale 
of the local instability. 

For finite times, absolute distinction between 
regular and chaotic orbits based on their Liapunov 
exponents is impossible. It is however the ex- 
ponential timescale that is important in discern- 
ing the physical significance of the instability. It 
is sufficient, for our purposes, for the inverse of 
the Liapunov exponent to be larger than a Hub- 
ble time for a trajectory to be considered stable. 
Elementary estimates and test calculations show 
that it requires about 50, 000 Myr for the value 
of the maximal exponent of trajectories to reach 
the value of about 10~ 4 Myr -1 . Suppose, for ex- 
ample, a flat rotation curve with a corresponding 
rotation period T = ^- — 60 Myr. In this case the 
divergence between neighboring circular orbits is 
SX ~ tot <~ 0. It and the exponents after 50,000 
Myr is - In 5, 000/50, 000 - 10~ 4 . 

The procedure of integrating a trajectory for 
far more than the age of the Universe to test 
whether it is stable on a smaller timescale can 
be rationalized in the following manner. For one 
can show that in a model barred potential most 
chaotic orbits can be separated from regular ones 
by inspection of their "local" Liapunov indicator 
a,i = In , where & is an initial infinitesimal 
perturbation and is its consequent on a sur- 
face of section, averaged over 10-20 consequents 
(Patsis et al. 1997). Note that except for the 
averaging over consequents and not over a fixed 
time interval, which has no effect on the result, 
the above scheme is equivalent to taking low n 
sums of Eq. (A12). In this context then, one can 
divide each of our trajectories into much smaller 
segments, say of 1 Gyr time-length each. Next, 
one can assume that each segment represent a new 
trajectory with initial conditions corresponding to 
the terminal end of the previous segment. Within 
this framework, the time-dependent Liapunov ex- 
ponent of this "mega-trajectory" can be consid- 
ered an ensemble average over many trajectories 
integrated over much smaller time. 

Because the exponential divergence can lead to 
mixing (as in a drop of ink mixing in a glass of 
water), and because phase space density is con- 
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served along trajectories, an initial localized prob- 
ability distribution will diffuse and eventually fill 
the whole connected region with equal probabil- 
ity (see, e.g., Lichtenberg & Lieberman 1995, and 
in the context of galactic dynamics Kandrup 1998 
and Merritt 1999). The rate at which this hap- 
pens is characteristic of the connected region and 
is determined by its Kolmogorov entropy — a mea- 
sure of the rate at which the coarse-grained phase 
space volume occupied by such a group of tra- 
jectory increases (see Appendix). When mixing 
is efficient one can replace long time averages of 
quantities over a single trajectory by short time 
averages over ensembles of trajectories. In general, 
because of the existence of semi-permeable phase 
space barriers and other complications (see, e.g., 
Wiggins 1991), such assertions are hard to prove 
in realistic dynamical systems and may be diffi- 
cult to ascertain numerically (because trajectories 
may take very long before they fill their allowed 
phase space region). Nevertheless, test results by 
Kandrup & Mahon (1994) suggest its basic plau- 
sibility by showing that the local Liapunov indica- 
tors mentioned above, averaged over ensembles of 
initial conditions in a connected phase space do- 
main, approximate well the long time Liapunov 
exponent of a single trajectory of the same region. 

The configuration space density distribution 
calculated over long timescales can also be thought 
to correspond to that of an ensemble of trajecto- 
ries integrated over a much shorter time. Because 
of the qualifications mentioned above however, it 
is still important to check the effects of the in- 
tegration timescale on the distribution. This is 
especially true when examining transient effects, 
like the time evolution of trajectories initially dis- 
tributed near a symmetry plane into a three di- 
mensional distribution. This will done in Sec- 
tion 6.3. 

Finally, a word about the calculation and in- 
terpretation of the full set of exponents. The ex- 
ponential instability, when present, will lead any 
perturbation to align itself with maximal direc- 
tion of expansion in phase space; the most un- 
stable direction. Thus calculation of Liapunov 
exponents from a random perturbation will re- 
sult in the maximal one being found. In order 
to find the other exponents, one starts from a set 
of orthogonal perturbations in phase space and re- 
orthogonalizc them at chosen intervals, to prevent 



them from re-aligning themselves again with the 
direction of maximal expansion. Several meth- 
ods have been proposed for this procedure (e.g., 
Eckmann & Ruclle 1985; Wolf et al. 1985). The 
method used here invokes the Gramm-Schmidt al- 
gorithm and is described in the appendix. For 
systems of three degrees of freedom six exponents 
are found. Since Hamiltonian systems have a skew 
symmetric structure 2 they come in pairs of posi- 
tive and negative ones, meaning that to each di- 
rection of expansion there is a corresponding con- 
traction (hence the conservation of Poincare's in- 
variants and phase space volume). It follows that 
zero exponents come in pairs and correspond to 
integrals of motion. In time-independent poten- 
tials there are always two zero exponents, corre- 
sponding to energy conservation. Other conserved 
quantities correspond to additional pairs of zero 
exponents. Trajectories with zero exponents are 
regular, they conserve three integrals of motion. 

3. Configuration space grid 

For the purpose of examining the configuration 
space distributions of our computed trajectories 
we have devised alOlxlOxll cylindrical grid of 
bins (IR, 19, IZ) defined as follows: 



IR = INT 



100 



log 10 (fl + l) 
logi (-R max + 1) 



+ 1 



(2) 



where R ma x is the maximum grid size taken to 
be double the bar's major axis and the last bin 
includes all radii lying beyond it, 



IZ = INT 



10 



lo; 



»g 10 \z\ x 10+ 1) 
\og 10 (z max x 10+ 1) 



(3) 



where z max is the maximum vertical grid size 
taken to be six times the vertical height of the bar, 
and again includes all the space above it. Since the 
problem is symmetric with respect to the x — y- 
planc, only positive values need to be recorded. 
Finally, 



Id = INT 



tan" 1 



5-7T 



+ 1 



(4) 



where x and y are Cartesian coordinates in a frame 
rotating with the bar. Again, by virtue of sym- 
metry with respect to the x and y axes only one 



2 Hamilton's equations are symmetric in the variables except 
for the minus sign in one of them 
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quadrant need be considered (hence the absolute 
value). 

This frame provides high resolution in the cen- 
ter without being too coarse in the outer regions. 
This is necessary, as we will be considering tra- 
jectories ranging from 60 pc to 6 kpc in radius, 
and starting from 20 pc in the vertical direction, 
but venturing 3-4 kpc out. The binning interval is 
taken to be 0.05 Myr. 

4. Units and mass models 

We use the kpc as length unit and the Myr as 
the time unit. This fixes the mass unit at about 
2.2 x 10 11 M©, when the gravitational constant is 
taken as unity. 

We follow the traditional route of superposing 
distinct functional forms for each of the compo- 
nents, constructing our galactic models from sep- 
arate disk, bar and halo contributions. This pro- 
cedure is an idealization, made possible by the 
linearity of the Poisson equation and the corre- 
sponding additivity of the forces due to different 
components. It allows one to distinguish the dif- 
ferent conceptual components. The equations of 
motion, however, are nonlinear and the effects on 
the motion due to the different components can- 
not, in general, be disentangled (since, obviously, 
the motion in a disk-halo system, for example, is 
not simply the motion in the disk potential super- 
posed onto the motion in the halo potential) . Our 
components are thus not to be thought as approx- 
imating the density distributions arising from dif- 
ferent kinematical components in a galaxy, which 
arise from the nonlinear motion, by a linear super- 
position of uncoupled auxiliary components (sub- 
stitutes). The goal is to mimic in a reasonable way 
the galactic mass distribution, especially those as- 
pects most important to our analysis — radial 
density profile and symmetry properties. 

A well known simple expression for a triaxial 
figure which has an asymptotically flat rotation 
curve is 

®h = \v 2 m log (R 2 h0 + x 2 + py 2 + qz 2 ) , (5) 

where Vho is the maximum rotation velocity (in- 
variably taken to be 0.2 ~ 200 km s _1 ) and Rho is 
the radius of a harmonic core. The parameters p 
and q are related to the (ellipsoidal) equipotential 



axis ratios by p = (ah/bh) 2 and q — (ah/ch) 2 , with 
a-h > bh > Ch for a triaxial halo. In all our runs we 
fix Ch/a,h = 0.8, while bh/ah is varied in the range 
between 0.9 and 1. 

The azimuthally-averaged density varies with 
radius as <~ R~ a with radius, with a approach- 
ing when R < R h0 and 2 for R > R h0 . It 
is intermediate at radii of order R h0 . At these 
radii the density approximates that found at the 
inner and intermediate regions of halos identified 
in cosmological simulations of the cold dark mat- 
ter structure formation scenario. Fig. 1 compares 
the density profile of the logarithmic potential to 
the closest (by inspection) models with inner den- 
sity falling as 1/R and 1/R 1 - 5 — corresponding to 
inner density profiles found by Navarro, Frenk & 
White 1997 and Moore et al. 1999, respectively. 
In both these cases, the density falls off as 1/R 3 
in the outer regions. The cosmological halos have 
scalelength, determining the transition from the 
inner profile to the outer one, of 3.7Rho and WRho 
respectively. 




LOG R 

Fig. 1. — Density of logarithmic potential (solid 
lines) as function of radius (in units of the core 
radius i?o) as compared to those of closest (by 
inspection) NFW model with central density in- 
creasing as 1/R (top), and more concentrated sys- 
tems with density increasing as 1/R 1 ' 5 (bottom). 

The main role of the disk component is to intro- 
duce strong asymmetry normal to its plane and to 
reduce it in that plane. Although an exponential 
disk is most realistic, we have decided to use the 
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Miyamoto-Nagai (1975) model, because in some 
cases we will be comparing our results with the 
self-consistent models of Pfenniger (1984b), which 
had this form for the disk potential. It is also 
easily handled computationally with the force ob- 
tained in closed form by taking the derivatives of 
the potential given by 



\Jx 2 + y 2 + (a d + y/bj+z^) 

where GMn is the mass of the disk, in units of G = 
1. The parameters ad and bd determine the central 
concentration and the flatness of the disks. They 
are taken to have values of 2.5 and 0.5 respectively. 
The disk mass is always taken to be GMjj = 0.2. 

A rapidly rotating bar with a pattern speed Qfc 
has been added in the plane containing the halo 
minor axis. This plane also contains the disk when 
it is present. The bar model used was a Ferrers 
ellipsoid of order two (Binney & Tremaine 1987; 
Pfenniger 1984a). The density of the Ferrers bar 
is constant along contours given by 

m 2 = x 2 /a 2 b +y 2 /b 2 + z 2 /c 2 , (7) 

where a& > &b > c& are the semi-axes of the density 
distribution. Inside the mass distribution (m < 1) 
the density varies as 

p = Pc (l-m 2 ) 2 , (8) 

where the central density p c is determined by the 
total mass of the bar GMb- The density is zero 
for m > 1 . 

The bar semi-axes are taken as cib = 6, b b = 1.5 
and Cb — 0.5, for comparison with Pfenniger's 
(1984a, b) main model. The bar pattern speed was 
always chosen so that the bar extends to corota- 
tion — calculated while assuming an axisymmetric 
halo and not including the bar's own contribution 
to the potential. Thus in practice the bar ends 
inside corotation in accordance with empirical re- 
quirement (e.g., Athanassoula 1992). In models 
where a disk was present, the bar had 20% of the 
disk mass within the circle tangent to the edge of 
its major axis. 

Models 1 — 3 include all the components, the 
triaxial halo, bar and disk. They are ordered in se- 
quence of increasing halo concentration. Model 1 



is a "maximal disk" model while Model 3 is halo 
dominated at all radii. The rotation curves of 
these three models are shown in Fig 2. Mod- 
els 4 — 6 arc without a disk. This makes it 
straightforward to quantify the net contribution of 
the non-rotating non-axisymmetric perturbation 
- the (constant) cquipotential axis ratios mea- 
sured in the inertial frame in the absence of the 
bar are simply those of the halo. Model 7 has 
an axisymmetric halo and a bar only, Model 8 
corresponds to main model studied by Pfenniger 
(1984a) with a bar embedded in a disk without 
halo. Model 9 deals with triaxial halo only. The 
complete set of model parameters is given in Ta- 
bic 1. 

5. Coordinate system, initial conditions 
and integrator 

We employ the coordinate system of Pfenniger 
(1984a), whereas the canonical phase space vari- 
ables are given by (x,y,z,X,Y,Z) and through these 
the Hamiltonian is defined as 

H = \ (X 2 + Y 2 + Z 2 )+<S>{x, y, z)-n b (xY - yX) , 

(9) 

where the potential $ is given by the sum of the 
contributions described in the previous section. 
In this system the spatial coordinates refer to a 
frame uniformly rotating with the bar with an- 
gular velocity fl bl while the velocities are mea- 
sured in the inertial frame. The coordinates of 
the halo, assumed to be non-rotating, are rotated 
back according to Xh — x cos(tt b t) — y sm(fl b t) and 
yh = xsm(Q,i,t) + y cos(Qfot). Once calculated, the 
halo force is then transformed back into the rotat- 
ing frame in a similar manner. 

Our goal is to investigate the stability of candi- 
date trajectories which may support a bar embed- 
ded in different disk-halo configurations. For this 
purpose the following initial conditions are appro- 
priate. We start slightly above the x-axis, directed 
along the bar major axis. The initial amplitude 
of the z-offset is taken as 20 pc for runs, except 
for those in Section 6.3 where the initial z excur- 
sion is increased to 200 pc. The velocities of the 
particles are taken normal to the xz-plane: that 
is the only nonzero component is Y, also taken 
to be positive. This means that our trajectories 
are symmetric with respect to the x-axis and they 
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Table 1 
Model Parameters 



Model 


RhO 


Vho 


bh/a h 


GM D 


ad 


b d 


GM B 




Notes 


1 


10.0 


0.2 


0.9 


0.2 


2.5 


0.5 


0.029 


0.031 


Triaxial halo with disk, bar 


2 


2.0 


0.2 


0.9 


0.2 


2.5 


0.5 


0.029 


0.041 


Triaxial halo with disk, bar 


3 


0.5 


0.2 


0.9 


0.2 


2.5 


0.5 


0.029 


0.042 
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are also, at least in the absence of additional tri- 
axial perturbation, prograde. In the case where 
the halo is axisymmetric, these initial conditions, 
within the limit of resolution in the initial x and 
Y values, produce all trajectories parented by the 
x\ periodic orbits aligned with the bar. These are 
known to be the building blocks of self-consistent 
bars (e.g., Pfenniger 1984b, where they are termed 
B orbits). 

The binning of the initial conditions is done on 
a linear equal-spacing grids with 100 subdivisions 
in each of these coordinates. Thus in total we 
have, for each model, 10,000 trajectories, with 100 
of them starting from each x position. The maxi- 
mum x position is taken to be the bar major axis 
(6 kpc) and the maximal Y velocity at each radius 
is 1.25 times the local rotation speed assuming an 
axisymmetric halo (and neglecting the bar's own 
contribution) . 

In some of the figures trajectories will be la- 
beled sequentially in the following manner. The 
trajectories with the lowest initial x coordinate are 
taken to be the first hundred. They are ordered in 
ascending manner according to their value of Y . 
Thus the trajectory having the lowest rotation ve- 
locity is number one and that having the highest is 
number 100. Next come the trajectories with the 
second lowest initial x values, again ranked in as- 
cending order, according to their initial Y values. 
And so on in a way that the last hundred (those 
with rank 9,900 to 10,000) trajectories start with 



largest initial x values in ascending order in Y. 

The integration is advanced using the variable 
order variable stepsize Adams method as imple- 
mented in the NAG routine D02CJF with a lo- 
cal (per timestep) tolerance of 10~ 10 . Even at 
this low tolerance level there is no guarantee that 
the chaotic trajectories integrated are the actual 
trajectories of the system from the given initial 
conditions. Indeed, varying the tolerance level 
gave usually different trajectories. This is be- 
cause the problem is inherently unstable. Nev- 
ertheless, it was found that for the range of toler- 
ance 10~ 8 — 10~ 14 , trajectories were qualitatively 
similar, as evidenced by their occupation numbers 
on the grid of Section 3, for example. In addi- 
tion, when canonical conjugate initial conditions 
were chosen for the tangent space vectors (i.e., vec- 
tors pointed towards pairs of position and velocity 
coordinates (x,X), etc.), the Liapunov exponents 
came in positive and negative pairs to an accu- 
racy of better than a few parts in 10 8 . This is a 
good measure as to the accuracy of the calcula- 
tion since it implies that the Poincare' invariants 
of all order (e.g., Arnold 1989; Sussman, Wisdom 
& Mayer 2001) are conserved to high accuracy and 
the symplectic nature of the Hamiltonian system 
is thus conserved: despite the fact that the inte- 
grator is not symplectic by design, it is effectively 
so. 

Except for Section 6.3, where results are pre- 
sented for trajectories integrated through a time 
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Fig. 2. — Rotation curves of Model 1 (top), 
Model 2 (center) and Model 3. Dashed curves rep- 
resent halo contributions, dashed dotted curves 
are the disk rotation curves while the solid lines 
relate to the total rotation curves. 

interval of 10, 000 Myr, all data relate to trajecto- 
ries evolved for 50, 000 Myr. Arguments rational- 
izing the choice of this particular time interval are 
given in Section 2. 

6. Stability of trajectories and configura- 
tion space distribution 

In this section we will analyze the stability 
properties of trajectories and their relationship 
to trajectories distribution in the configuration 



space. This will be achieved by calculating the 
Liapunov exponents and by employing the cylin- 
drical grid described in Section 3. Our goal will be 
to examine the plausibility of self-consistent equi- 
libria for barred galaxy models with different halo 
structures — differing from each other in their de- 
gree of central concentrations and triaxiality. In 
addition, we would like to know how tight is the 
correlation between the trajectories' distribution 
in the configuration space and their stability prop- 
erties. 

6.1. Maximal exponent and configuration 
space volume 

6.1.1. The effect of central concentration 

To get a measure of the configuration space vol- 
ume occupied by a given trajectory we calculate 
the total volume of the cells which it visits. This 
volume is normalized in terms of the bar's volume 
within the trajectory starting point, i.e., the vol- 
ume of the bar between two planes normal to the 
major axis of the bar at the starting point (and its 
reflection). The results for Models 1 — 3 are shown 
on the left-hand panels in Fig. 3, which depict 
a sequence of increasingly centrally-concentrated 
halos with core radii decreasing from 10 kpc to 
0.5 kpc. The models are presented in the form 
of grayshade diagrams, with darker shading corre- 
sponding to larger relative volume occupied by a 
trajectory. The positions on the diagrams corre- 
spond to the initial conditions — with the vertical 
axis rescaled in terms of the local rotational speed 
at the initial x value labeled on the horizontal axis. 

The darkest shading corresponds to trajectories 
occupying a configuration space volume which is 
100 or more in the normalized units (see above). 
Clearly such trajectories (and all those that are 
represented by the darker spots on this logarithmic 
scale) cannot contribute towards a self-consistent 
bar — by virtue of the fact that their spatial dis- 
tribution is far more extended than the bar at 
their starting radius. As can be seen, the fraction 
of such trajectories increases significantly as the 
halo central concentration is increased: being con- 
fined to the outer regions (where the equidensity 
profiles of the model are themselves more round) 
in Model 1, but occupying the vast majority of 
all initial conditions in Model 3. As we will see 
shortly, these trajectories are chaotic, with large 
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Fig. 3. — Grayshade diagrams showing configura- 
tion space volume (left panels) and maximal Lia- 
punov exponents (right panels) for trajectories of 
Model 1 (top), Model 2 (middle) and Model 3. 
The configuration space volume is defined as the 
total volume of the cells visited by a trajectory 
and is normalized in units of the volume of the 
bar within the starting radius of the trajectory. 
The scale is logarithmic (base ten) with the white 
background corresponding to values of —0.5 and 
less, and the foreground to values of 2 (i.e., hun- 
dred times the volume of the bar at the initial 
radius) and more. The exponents are rescaled, so 
that their inverse is given in units of dynamical 
time (taken to be the local rotation period in the 
azimuthally-averaged potential excluding the bar) . 
The scale is logarithmic. The background corre- 
sponds to values of —1.5 or less (i.e., exponential 
times of ~ 32t£> or more and the foreground to 
values of 0.5 or more. 

Liapunov exponents. And even though the vol- 
ume of configuration space occupied by some of 
the chaotic trajectories in this model is sometimes 
smaller than corresponding trajectories in the case 



of Rha — 2, trajectories with volume significantly 
larger than the bar are much more abundant. In- 
deed, in this latter case, except for the "island" 
around V y /V ro t = 0.75 between initial condition 
with 1.5 kpc and 3 kpc and the very inner region 
(where the halo potential is nearly harmonic) there 
are hardly any regular trajectories that are elon- 
gated with the bar Moreover, many of the regular 
trajectories do not match the bar density distribu- 
tion, being too round by comparison. This is the 
case of the regular regions corresponding to initial 
velocities Vy/V rot > 0.8 and initial radii greater 
than 3 kpc. For as can be seen from Fig. 4, where 
we plot the distribution of the minimum (cylindri- 
cal) radius visited by these trajectories, the major- 
ity of trajectories started beyond x = 3 have min- 
imal radii that are greater than the bars minor 
axis; they envelop the bar and thus cannot con- 
tribute towards building a self-consistent model. 
Note that the trajectories which do have minimal 
radii smaller than 1.5 correspond to the isolated 
islands of stability corresponding to initial veloc- 
ities Vy/V ro t iS 0.8. One can, therefore, conclude 
that the sequence of Models 1 — 3 represents pro- 
gressively more unstable bars which will tend to 
evolve quickly, since the system does not support 
the kind of trajectories required to build such a 
bar self-consistcntly. 

From the right hand-side panels of Fig. 3 one 
observes that most of the trajectories occupy- 
ing large regions of the configuration space are 
chaotic. Furthermore, comparison of left and right 
panels shows that there is quite a tight correla- 
tion between the values of the Liapunov exponents 
and the volume of configuration space the trajec- 
tories move in. This, of course, should not come 
as a surprise. For as explained in the introduction 
and in Section 2, regular orbits are confined by 
invariants which characterize the regular motion 
while chaotic ones are under no such constraints. 
Trajectories with higher Liapunov exponents, in 
addition to having smaller instability timescales, 
are also typically found in regions of phase space 
away from regular trajectories and barriers, such 
as cantori which can prevent phase space transport 
(e.g., Wiggins 1991). This implies that, in a given 
time, they are to diffuse in larger regions of phase 
space, which is then reflected in their configura- 
tions space distribution. Nevertheless, the degree 
of correlation is stark and shows the intimate rela- 
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Fig. 4. — Distribution of minimal radii of regu- 
lar trajectories (defined here as those having con- 
figuration volume smaller than the bar's at their 
starting radius) of Model 3 with initial coordinate 
x > 3. 

tion the question of stability of trajectories bears 
to that of self-consistency of galactic models. 

There are two effects which may be leading to 
the increase in the fraction of chaotic trajectories 
with the decrease in the halo core radius. The 
first, as discussed in the introduction, is the in- 
creased nonlincarity introduced by the centrally- 
concentrated mass distribution — which leads to 
increased coupling between the degrees of freedom. 
The second is the increased time dependency in 
the force field. For, as can be seen from Fig. 5, 
where we plot the ratios of the non-integrable com- 
ponents of the bar and halo force fields, for the 
least centrally-concentrated halo (Rho = 10) it 
is only at the edge of the bar that the average 
azimuthal force component of the halo becomes 
comparable to that of the bar. In the case of 
Rho = 0.5, on the other hand, this component is of 
the order of the corresponding bar component at 
all radii. For mild halo triaxialities, therefore, it is 
necessary that the halo be centrally-concentrated 
for the effect of time dependency to become im- 
portant. 




RADIUS 

Fig. 5. — Ratio of the average of the absolute value 
of the bar's azimuthal force to that of the halo, 
for halo axis ratio b/a = 0.9 and for different core 
radii. The bar's relative contribution is seen to 
be strongest when the halo has a large core radius 
i?o = 10, weakest for a centrally-concentrated halo 
R = 0.5 and intermediate for R = 2. 

6.1.2. The effect of halo triaxiality 

We now attempt to look in more detail into 
the origin of the widespread chaotic behavior ob- 
served in the orbital structure of barred systems 
with centrally-concentrated halo. To more eas- 
ily identify the strength of the non-rotating non- 
axisymmetric contribution to the potential, we 
now remove the axisymmetric disk from our su- 
perposition of potentials. Removing the disk con- 
tribution will obviously increase the effective non- 
axisymmetric perturbation in the potential. This 
is expected to increase the region of instability 
even further. It is apparent, however, that the 
dominant element producing unstable trajectories 
near the short axis plane of centrally-concentrated 
halos in our models is not their own triaxiality 
but the presence of a rotating barred component. 
This can be seen from Fig. 6, where we show the 
grayshade diagrams for a centrally-concentrated 
triaxial halo with b^ja^ = 0.9, Ch/dh = 0.8. In 
the absence of the bar, trajectories are regular and 
are confined in configuration space (Model 9, top 
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panels). The addition of the rotating perturba- 
tion, however, modifies the situation dramatically 
(Model 4, center). For comparison, we show the 
grayshade diagrams of the model self-consistcntly 
constructed by Pfenniger (1984b: Model 8). It is 
evident that the orbital structure, as manifested 
by the configuration space distribution and Lia- 
punov numbers, of this model is different from 
that of models where a centrally-concentrated tri- 
axial halo is present. This of course confirms 
our original inference that a self-consistent barred 
system is impossible to sustain inside a triaxial 
and centrally-concentrated halo: the trajectories 
it supports are too chaotic and occupy too large 
a volume of the configuration space to represent a 
bar. 




Fig. 6. — Same as in Fig. 3 but for trajectories of 
Model 9 (top), Model 4 (middle) and Model 8. 

In fact, as will be seen from Fig. 7 (Model 6, 
middle panels), even very small departures from 
axisymmetry on the part of the halo (as small as 
1% in the potential axis ratio bh/ah) can lead to 
large connected regions of initial conditions corre- 
sponding to chaotic trajectories, with only islands 



of stability in their midst. Only in the limit of 
a perfectly axisymmetric halo models is it possi- 
ble to obtain connected regions of regular trajec- 
tories at most radii. In that case (bottom pan- 
els), significant fraction of the initial conditions 
are occupied by trajectories that support the bar 
— the white regions that appear for most initial x 
up to about 4 kpc represent trajectories that are 
mostly parented by stable periodic x\ orbits, and 
are thus aligned with the bar. 3 The existence 
of these trajectories is necessary for the construc- 
tion of barred equilibria. It is consistent with nu- 
merical simulations where long lived bars embed- 
ded in centrally-concentrated axisymmetric halos 
are observed (Ideta & Hozumi 2000; Debattista & 
Sellwood 2000; Athanassoula & Misiriotis 2002). 
These may slowdown in time but do not dissolve. 
In some of these simulations rings surrounding the 
bars are also found. These are possibly related to 
the trajectories represented by the white strips on 
the upper right corner of the diagrams in the bot- 
tom panels of Fig. 7. These regions correspond 
to nearly round trajectories which, at the outer 
edge of the bar and slightly beyond, are elongated 
in the direction perpendicular to the bar — as is 
observed in the simulations. 

The fraction of bar supporting regular trajecto- 
ries decreases rapidly as small non-axisymmetric 
perturbations are applied. In the case of Model 6 
orbits elongated along the bar still exist at most 
radii. However they occupy a much smaller frac- 
tion of the initial conditions than in Model 7. The 
range of shapes they come in is therefore narrower, 
which implies that it is less probable that among 
them will be found those that match the bar's 
asymmetry (i.e., are not too fat). Indeed, we find 
by inspection that a large fraction of the avail- 
able trajectories that are elongated along the bar 
do not match its asymmetry. Further departures 
from axisymmetry render the bar supporting or- 
bits exceedingly rare. In Model 5 (Fig. 7, top pan- 
els) they are found at a small starting range in x. 
In model 4, they are virtually non-existent. 

The above implies that our axisymmetric model 

3 The periodic x\ orbits survive, in distorted but stable form, 
for mild halo triaxialitics (b h /a h > 0.92). For potential axis 
ratios of 0.9 or smaller, they are (at least mostly) unsta- 
ble. Detailed examination of the existence and stability of 
these periodic orbits in time-dependent potentials will be 
presented elsewhere. 
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Fig. 7. — Same as in Fig. 3 but for trajectories of 
Model 5 (top), Model 6 (middle) and Model 7. 

is structurally unstable — since any small pertur- 
bation renders the once connected regular regions 
into disjoint sets. The replacement of regular tra- 
jectories parented by the X\ orbits by chaotic ones 
in turn implies that self-consistent bars are un- 
likely under these circumstances. Because the 
perturbations required are very small, axisym- 
metric models of barred galaxies with centrally- 
concentrated halos must then be considered non- 
generic. 

6.2. Vertical stability of trajectories 

The conclusion that systems with centrally- 
concentrated triaxial halos are unlikely supporters 
of bars is impressed further when one looks at the 
vertical stability of the trajectories under consid- 
eration. In Fig. 8 we show, for some of our mod- 
els, grayshade diagrams representing the fraction 
of time trajectories spend in cells that correspond 
to vertical excursions greater than the semi-minor 



of the bar. It is evident that for models with even 
minor triaxiality this fraction can exceed one-half 
of the time for a large number of trajectories. Note 
that this is necessarily an underestimate — the 
bar's extension is, except at the origin, always less 
than the minor axis value. 

Comparison of the structure of the grayshade 
diagrams in Fig. 8 with those in Fig. 6 and Fig. 7 
reveals that the fraction of time a trajectory 
spends having large absolute values for its ver- 
tical coordinate correlates well with the values of 
the Liapunov exponents. In other words, trajec- 
tories which spend the largest fraction of time at 
large z have largest Liapunov exponents. This 
in turn correlated with the total configuration 
volume they occupy. Thus vertical stability is a 
simple diagnostic of chaotic behavior: almost all 
regular trajectories conserve their initially small 
amplitudes of z oscillations (the exceptions are 
those near bifurcations leading to three dimen- 
sional regular orbits). It is also a sufficient crite- 
rion that orbit densities do not match that of the 
bar, being too extended in the vertical direction. 

We find that the maximal vertical excursion of 
chaotic trajectories of some of our models can be 
quite large — of the order of 10 kpc. To illustrate 
this, we have arranged the trajectories in ascend- 
ing order according to their starting spatial and 
velocity coordinates (see Section 2) and plotted 
the logarithm of their maximal vertical excursion 
within the first 10, 000 Myr. The results are shown 
in Fig. 9. 

6.3. Evolution on shorter timescale and 
the instability in- and out-of-the 
plane 

Some trajectories which have large Liapunov 
exponents, and according to the diagrams of Fig. 8 
spend a large fraction of their time at z coordi- 
nates larger than the bar's vertical extension, nev- 
ertheless appear to have small maximal vertical 
excursion in Fig. 9. This happens because many 
of these trajectories do not, for the small initial 
z amplitude we used, reach their maximal z ex- 
tension in 10, 000 Myr. We had chosen the small 
initial z amplitudes in order to detect trajecto- 
ries that may be vertically unstable even if their 
maximal excursions do not exceed a few hundred 
pc. In realistic situations however, vertical mo- 
tions would be sufficient so that stars already have 
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Fig. 8. — Grayshade diagrams showing the frac- 
tion of time trajectories spend at cells compris- 
ing vertical coordinates larger than the bar minor 
axis. Left panel: Model 9 (top), Model 4 (middle) 
and Model 8 (as in panels of Fig 6). Right panel: 
Model 5 (top), Model 6 (middle) and Model 7 (as 
in panels of Fig 7). The shades are on a linear 
scale whereas the foreground corresponds to tra- 
jectories spending 75% of their time or longer in 
the aforementioned cells and the background cor- 
responds to trajectories spending 0.075% of the 
integration time in such cells. 

z amplitudes of such magnitude. 

We have rerun some models with the larger ini- 
tial z amplitude of 0.2. The results, for Models 4, 5 
and 6 are shown in Fig. 10. One finds that many 
trajectories, that apparently were stable in the z 
direction over the Hubble time, arc now unsta- 
ble. Fig. 10 (right panels) also displays diagrams, 
similar to those in Fig. 8, representing the frac- 
tion of time trajectories spend in cells comprising 
vertical coordinates larger than that of the bar 



Fig. 9. — Maximal vertical excursion within 
10, 000 Myr. The panels correspond to those of 
models shown in Fig. 8. The trajectories are 
arranged in ascending order such that those at 
smaller initial radii have strictly smaller number 
than those starting at a larger radius. At each ra- 
dius the rank progressively increases with initial 
normal velocity (cf., Section 5). 

vertical extension. This fraction of time now has 
increased significantly by starting our trajectories 
at a larger z coordinate value, albeit one that still 
is significantly smaller than the bar's semi-minor 
axis. This suggests that the vertical instability 
is important in determining the orbital properties 
of the system, transforming trajectories that are 
stable in the plane into unstable ones, but that 
it manifests itself over relatively long time-scale, 
especially when the initial perturbation is small. 

Trajectories having large vertical extensions 
and spending most of their time with absolute 
values of their z coordinates larger than the bar's 
semi-minor axis cannot contribute towards a self 
consistent barred galaxy model. This being the 
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Fig. 10. — Maximum vertical excursions (left pan- 
els) and grayshade diagrams showing the time 
fraction spent at configuration space cells com- 
prising vertical coordinates larger than the bar 
minor axis for Model 4 (top), Model 5 (middle) 
and Model 6 (with the same scaling as in Fig. 8. 
Trajectories are integrated for 10, 000 Myr start- 
ing with a vertical amplitude of 0.2 kpc (instead 
of the canonical 0.02 kpc used in the other runs). 

case for most trajectories of Models 4,5 and 6, 
makes it apparent that the bar would not sur- 
vive in such models for a Hubble time — which 
places an upper limit on the dissolution time of 
any barred system embedded in mildly triaxial 
centrally concentrated halo. A lower limit, which 
should be of the order of the Liapunov timescale 
of the chaotic trajectories, is of the order of a few 
dynamical times. This should be the time-scale 
for bar dissolution in the most chaotic systems 
(e.g., Model 4) where hardly any regular trajecto- 
ries exist. In such cases, the vertical instability is 
unlikely to be of central importance in the evolu- 
tion of the system — since the bar would probably 



dissolve in the plane before its presence becomes 
felt. 

Fig. 11 shows the maximal Liapunov exponents 
and the configuration space volume for Models 5, 
6 and 7 for trajectories integrated through this 
smaller period of time and starting with larger z 
amplitude (note the difference in grayscales of Lia- 
punov exponents). Again one observes a clear sim- 
ilarity between the distribution of values of the ex- 
ponents and the occupied configuration space vol- 
ume. The figures are qualitatively similar to the 
ones integrated for 50, 000 Myr from z in = 0.02 
(cf. Fig. 7). Quantitatively there are some dif- 
ferences. The extent of the shaded regions has in- 
creased. Unstable trajectories are now more abun- 
dant. This is again an effect of the larger ini- 
tial vertical amplitudes used here, causing chaotic 
trajectories to replace regular ones. However the 
most unstable trajectories now have somewhat 
smaller configuration space volume than previ- 
ously. This suggests that trajectories continue to 
explore new regions of phase space as the system 
evolves in time. It is not surprising, given the 
time-dependent nature of the potential (meaning 
an orbits phase available space is not necessarily 
bounded) and the possibility of the presence of 
"sticky" chaotic trajectories that take long times 
to achieve an invariant distribution even for the 
time- independent case (e.g., Siopis & Kandrup 
2002). 

Indeed simple inspection of the spatial distribu- 
tion of trajectories showed that both these effects 
are at work; some trajectories keep on moving to- 
wards larger and larger radii and vertical excur- 
sions, while others stay within a bounded region 
but fill larger and larger volumes within that re- 
gion. In the latter case this mainly consisted of 
trajectories with a hollow configuration space dis- 
tribution, which was being filled up as the inte- 
gration proceeds. There are also some trajectories 
whose vertical extension is only significant after 10 
Gyr. None of the trajectories examined however 
had a spatial configuration that did not match the 
bar over 50 Gyr but did match the bar over 10 Gyr. 
In fact, the vast majority already showed a spa- 
tial structure in the plane which is very different 
from that of the bar within 500 Myr or so — their 
shapes being too isotropic, even though they were 
launched from initial conditions corresponding to 
thin bar-aligned orbits in the case when the halo is 
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Fig. 11. — Same as in Fig. 7 but for trajectories 
integrated for 10, 000 Myr and with initial verti- 
cal condition z = 0.2 kpc (instead of the canonical 
0.02 kpc used in the runs in that figure). Note 
that, to adjust for the significantly shorter inte- 
gration time, which implies that regular trajecto- 
ries would have an exponent that is proportionally 
larger, the background for the grayshades repre- 
senting these quantities is rescaled to -0.9 (instead 
of -1.5 as in Fig. 7). 

axisymmetric. Vertical instability generally took 
longer to develop, but was far faster, as expected, 
when the initial vertical perturbation was larger. 

Furthermore, trajectories that occupied a rela- 
tively (compared to the bar) large configuration 
space over the longer times integration still do 
over the span of 10 Gyr. This can be seen from 
Fig. 12 where we plot the relative volumes of the 
trajectories of Model 6 for the two cases against 
the volumes of the case where the trajectories are 
integrated for 50 Gyr. As expected, trajectories 
which, for the case of long integration from small 
initial z, had small phase space volumes usually 



acquire a much larger one. This is due to the 
enhanced z instability which further destroys the 
regular regions. On the other hand, trajectories 
which had large configuration space volumes com- 
pared to that of the bar and had pronounced z 
excursion, even when starting from z = 0.02, now 
have smaller volume. Nevertheless, this is still 
much larger than the that of the bar — again 
confirming that they are unlikely to contribute to- 
wards a self-consistent model on the timescale of 
10, 000 Myr. 




Fig. 12. — Comparison of the configuration space 
volume for trajectories integrated for 10000 Myr 
(and stating AT z = 0.2) and ones integrated 
for 50000 Myr (and starting at z = 0.02). The 
the horizontal axis represent values of the (base 
ten) logarithm of the configuration space volume 
of the latter trajectories, while the vertical axis 
represents logarithms (again base ten) of the ratio 
of the two volumes for the different trajectories 
of Model 5 (i.e., wo/(10000, i)/wo/(50000, i), i = 
1,10000). 

6.4. Full set of Liapunov exponents and 
distribution of exponential times 

Because it is CPU-time consuming to calculate 
the full set of Liapunov exponents at high reso- 
lution in the initial conditions, we performed this 
task for a selected set of models. The values of 
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these exponents provide information on whether 
trajectories, though being chaotic, may conserve 
(even if approximately) some quantities with time. 
Fig. 13 shows grayshade diagrams involving the 
second and third Liapunov exponents of trajecto- 
ries in the three models, 5, 6 and 7, whose config- 
uration space volume and maximal exponents arc 
shown in Fig. 7. The agreement between the con- 
figuration volume of trajectories and the values of 
their exponents is even more precise than in the 
case of the maximal exponent — most discrepan- 
cies, even though minor originally, having now dis- 
appeared. The additional constraints provided by 
low values in the two smaller Liapunov exponents 
account for the remaining regions of initial con- 
ditions corresponding to small configuration space 
volume despite having large values of the maximal 
exponent. 

The sum of the positive exponents for a given 
trajectory characterizes its Kolmogorov entropy 
(see Appendix). For a group of nearby trajecto- 
ries in connected phase space region it measures 
the rate of increase of the coarse-grained volume 
they occupy, and thus of the evolution of their sta- 
tistical distribution. The inverse of this entropy 
(strictly speaking multiplied by a numerical fac- 
tor of order unity that is ignored here) defines a 
timescale for such evolution. 

In Fig. 14 we display histograms showing the 
distribution of the inverse of the Kolmogorov en- 
tropy for trajectories of Models 7, 6, 5 and 4. 
As can be seen, for trajectories of axisymmet- 
ric Model 7, the corresponding exponential times 
are mostly very large, of the order of the Hub- 
ble time. These correspond mainly to trajecto- 
ries trapped around the x\ periodic orbit family 
aligned with the bar (the large connected white 
patch in the grayshade diagrams in the bottom 
panel of Fig. 7). Moreover, the more unstable tra- 
jectories, with smaller exponential timescales, are 
of roughly equal numbers over a large range of val- 
ues; the distribution is relatively flat, being only 
slightly bimodal on far left of the figure. These 
are exponents of various higher order regular or 
trapped chaotic trajectories — i.e., ones existing in 
regions of phase space where mainly regular orbits 
dominate. The final peak corresponds to trajec- 
tories starting from initial conditions in the outer 
region of the system (cf., Fig. 7 bottom panel), 
which correspond to regions of phase space that 
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Fig. 13. — Grayshade diagrams for the middle 
(left panels) and smallest Liapunov exponents for 
Model 5 (top), Model 6 (middle) and Model 7. On 
the (base ten) logarithmic scale the middle expo- 
nent background value is —1.5 and foreground is 
0.2 while for the smallest exponent the background 
is scaled to —2 and foreground to 0. 

are dominated by chaotic trajectories. 

As the non-axisymmetric perturbation is in- 
creased, regular trajectories aligned with the bar 
progressively become less abundant. Initially 
there is an equal increase in the number of trapped 
trajectories (second peak on left for histogram of 
Model 6) and highly chaotic ones (first peak). 
Eventually, however, the increase in the latter 
occurs at the expense of the former (Model 5). 
Finally, the fraction of both the trapped and reg- 
ular trajectories becomes very small (Model 4). 
Here most trajectories are highly chaotic, part of 
the same connected "stochastic sea," with expo- 
nentiation timescale varying relatively little. This 
is expected since, for Model 4 (cf. Fig. 6, middle 
panel), except for small regular regions, trajecto- 
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Fig. 14. — Histograms showing the distribution of 
the exponentiation times (taken to be the inverse 
of the Kolmogorov-Sinai entropy) for trajectories 
of (from top to bottom) models 7, 6, 5 and 4. 

ries form a connected region of initial conditions 
that have similar Liapunov exponents and config- 
uration space volumes. 

For Model 4, most of the remaining regular tra- 
jectories are parented by surviving higher order 
periodic orbit families of the resonant zones (the 
white strips in Fig. 6 middle panels). None are 
parented by the x\ family of the unperturbed ax- 
isymmetric halo case. Almost certainly, the ab- 
sence of a significant fraction of regular orbits, 
more importantly the total absence of any regular 
orbits aligned with the bar, ensures that such a bar 
cannot be built self-consistcntly. For in the case of 
a time-independent potential, for example, any set 
of initial conditions starting in a connected chaotic 
region will evolve to an invariant distribution, with 
every trajectory having the same time-averaged 



phase space distribution in as all other trajecto- 
ries. The projection onto configuration space of 
such a distribution would be too isotropic to match 
a bar. In the case of a time-dependent system, the 
region in which trajectories move is not necessar- 
ily bounded by a zero velocity curve, and therefore 
some trajectories can explore even larger volumes 
as their energies change — an effect that is likely 
to make the discrepancy more pronounced. 

7. Energy decorrelation 

For trajectories with large maximal exponent, 
the other two positive exponents are proportional 
to it with a nearly constant of proportionality. 
This is not always true, however, of trajectories 
with smaller maximal exponent. The latter tend 
to have relatively even smaller corresponding val- 
ues for the two smaller exponents. This is illus- 
trated in Fig 15 where we plot the Liapunov ex- 
ponents of trajectories of Model 4. It is especially 
clear for trajectories starting in the inner regions 
(those with small trajectory numbers). It turns 
out that such trajectories also conserve their en- 
ergy well. 

Energy, always an invariant quantity in a time- 
independent potential, is allowed to vary in the 
time-dependent models described here. Corre- 
sponding to energy conservation is a pair of zero 
Liapunov exponents. Therefore, trajectories that 
conserve energy well will have at least one small 
exponent. In order to test energy (i.e., the Hamil- 
tonian of Eq. [9]) conservation, we calculated the 
correlator C; = J^t^SVilf 1 w ■> where r is a time 
interval, taken as 500 Myr, and i is an integer. The 
quantity d as defined above is zero if there is com- 
plete loss of memory in energy over a timescale ri, 
and one if there is no loss whatsoever, — for ex- 
ample, if energy is conserved along a trajectory. It 
can also be seen as an angle between a normalized 
unit vector, consisting of the values of the energy 
at intervals ri, and another "delayed" vector with 
corresponding points delayed from the first vector 
also by an interval ri. 

In Fig. 16 we plot log 10 (l — Cj) for i = 1, 10, as 
well as the relative dispersion in energy along our 
enumerated trajectories of Model 4. As expected, 
in the region where two of the Liapunov exponents 
have very small values, energy correlation is large 
(that is near 1) and the dispersion is small. What 
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Fig. 15. — Liapunov exponents of Model 5 for tra- 
jectories arranged sequentially such that the first 
hundred correspond to initial conditions at the 
first position bin and progressively larger normal 
velocities (cf., Section 5). 

is somewhat surprising is that, for larger radii, the 
energies of the trajectories at different times ap- 
pear to be again correlated — this despite the fact 
that these trajectories belong to large connected 
regions with high values of the exponents (cf. Sec- 
tion 6.1.2). 

The variation of the dispersion and the quan- 
tities log 10 (l — Cj) as a function of radius appear 
to follow closely the strength of the bar pertur- 
bation (Fig. 17). This trend does not seem to be 
substantially altered by changing the asymmetry 
of the halo (we have checked this by comparing 
the results presented here with analogous ones for 
Models 5 and 6). Thus, for any non-axisymmetric 
perturbation between one and ten percent in the 
nonrotating potential, energy conservation of tra- 
jectories starting at radii where the bar asymmetry 
is maximal is affected considerably. For other tra- 
jectories, energy conservation is not significantly 
affected by such small perturbations. One can 
perhaps say that the value of their Hamiltonian 
is stable along their trajectories. 



Fig. 16. — Energy correlators d = 1 — 
<E(t) ] ><E(I+T)> av eraged over 50,000 Myr with 
interval ri and r = 500 MYr (Top: i = 1, Middle: 
i = 10), and relative energy dispersions (bottom) 
for Model 5. Trajectories are arranged sequen- 
tially such that the first hundred correspond to 
initial conditions at the first position bin and pro- 
gressively larger normal velocities (cf., Section 5). 



It is interesting that a time-dependent pertur- 
bation can have such a large effect on the stabil- 
ity properties of some trajectories and the manner 
in which their configuration space distribution be- 
haves, but causes very little change in energy. It 
illustrates the fragility of orbital structure of the 
unperturbed axisymmetric system. 

8. Conclusions 

In this paper we have analyzed the stability and 
configuration space properties of trajectories start- 
ing near the plane containing a bar in galaxies with 
different dark matter halo distributions. In partic- 
ular, we examine the variation of these properties 
for different concentrations (Section 6.1.1) and tri- 
axialities (Section 6.1.2). This is done by calculat- 
ing the Liapunov exponents, looking at the time- 
averaged spatial density distribution (Section 6) 
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Fig. 17. — Average (of absolute value of) az- 
imuthal bar force as a function of radius. 

and the vertical stability properties (Sections 6.2) 
of samples of 10, 000 trajectories starting near the 
plane containing the bar in each model. As a ref- 
erence we have compared these properties with 
those of the relatively stable model that has been 
built self-consistently by Pfenniger (1984a, b). For 
some of the models we also calculate the full set of 
Liapunov exponents and look at the distribution 
of exponential timescales inferred from the Kol- 
mogorov entropy (Section 6.4). In general, there 
is remarkable correlation between the local stabil- 
ity properties, as characterized by the exponents, 
and the the spatial distribution of trajectories, as 
quantified, for example, by the configuration space 
volume they fill — with the more unstable trajec- 
tories occupying much larger volume than more 
chaotic ones. 

For models with centrally-concentrated ha- 
los, we find a large fraction of trajectories to 
be chaotic. These fill a large (compared to the 
bar's) volume of the configuration space and spend 
most of the time at vertical extension greater 
than the maximum thickness of the bar. They 
are, therefore, unlikely to contribute towards a 
self-consistent model. In addition many of the 
trajectories that remain regular have minimum 
extensions in the plane larger than the bar's mid- 
dle axis. They also cannot contribute towards a 



self-consistent model (Section 6.1.1). 

The above is even true to some extent in the 
case when the halo is centrally-concentrated but 
axisymmetric. Still, in this case, there remains 
large connected sets of initial conditions corre- 
sponding to regular orbits which are aligned with 
the bar and could contribute to a self-consistent 
model. When, however, the halo is triaxial, 
chaotic trajectories become a dominant majority 
— even very small halo triaxiality (a 1% perturba- 
tion in the potential axis ratio) causes the chaotic 
regions to become connected and regular phase 
space regions to become islands of stability (Sec- 
tion 6.1.2). Such models are unlikely candidates 
for steady state self-consistent solutions. 

There are several reasons why centrally-concentrated 
halos, especially triaxial ones, contribute to the 
chaotic behavior displayed from most initial con- 
ditions. A centrally-concentrated mass distribu- 
tion has solutions for the Poisson equation that 
is far from quadratic, in the sense that an expan- 
sion of the potential power series has large terms 
beyond the quadratic. This produces a coupling 
between different degrees of freedom in equations 
of motions which are highly nonlinear. In the 
axisymmetric halo case these systems generically 
have no global integrals of motion other than en- 
ergy (in the rotating frame) and, being far from 
linear, arc, therefore, candidates for supporting 
large (phase space) regions of chaotic trajectories. 
When a nonrotating triaxial halo is added even 
energy is lost as an integral of motion. Since time 
translation symmetry is broken, such a system 
has no obvious symmetries and is likely to exhibit 
chaotic behavior from most initial conditions, as 
indeed is found to be the case. Moreover, in the 
case of a halo with a significant core, and for mild 
triaxiality, the non-axisymmetric halo perturba- 
tion is completely negligible in the bar region. 
This is another reason why central concentration 
triggers chaos; for centrally-concentrated halos of 
mild triaxiality the bar and halo azimuthal forces 
in the equatorial plane become comparable. There 
are no frames of reference where the system may 
be considered, to some approximation, stationary. 

It is significant that very small deviations from 
axisymmctry in the halo potential can produce 
drastic changes in phase space structure. Sys- 
tems where neither KAM stability (a property of 
separable systems: see, e.g., Arnold in Mackay & 
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Meiss 1987) guarantees that small perturbations 
do not drastically change phase space structure or 
where structural stability (a property of strongly 
chaotic systems, e.g., C systems: Anosov 1969) 
ensures the robustness of qualitative properties of 
trajectories, are candidates for such behavior. In 
other words, systems with a mixed phase space, 
where regular and chaotic trajectories coexist, are 
always liable for strong modification of qualitative 
properties by means of small perturbations. This 
appears to be the case for models of centrally- 
concentrated barred galaxies with axisymmetric 
halos. A small perturbation in the halo axis ra- 
tio can, as we have seen, have a significant effect 
(cf., Fig 7). These systems are said to be struc- 
turally unstable. Modeling barred galaxies with 
centrally concentrated halos while assuming that 
these halos are axisymmetric may, therefore, pro- 
duce non-generic results. 

Even though we find some correlation between 
the energy change of trajectories of triaxial sys- 
tems and the values of Liapunov exponents, many 
trajectories that are unstable over a dynamical 
time conserve energy quite well (Section 7). This 
was inferred by calculating the energy correlation 
over times up to 5 Gyr and the dispersion in its 
values averaged over 50 Gyr. We find that only 
for a small fraction of the chaotic trajectories does 
energy decorrelate over this period and the disper- 
sion in its values becomes large. This shows that 
long timescales for energy relaxation do not imply 
that evolution of other quantities does not hap- 
pen at a much larger rate, as is the case with the 
z instability, for example. 

Halos identified in cosmological simulations 
with CDM initial conditions are found to be both 
centrally-concentrated and triaxial. The results 
presented in this paper, therefore, suggest that 
the existence of such structures around galax- 
ies containing strong large-scale bars is probably 
ruled out. Furthermore, the very small departures 
required from axisymmetry suggest that our re- 
sults are generic — in the sense that departures 
from perfect axisymmetry, due to halo substruc- 
ture, for example (another prediction of CDM 
models), may be sufficient for the effects found in 
this work to become important. These findings 
can, therefore, be considered consistent with the 
growing observational evidence against such halos 
in present day galaxies. 



It has been suggested that primordial bars in 
centrally-concentrated halos can act as to reduce 
the central concentration of the halo (e.g., Binncy, 
Gerhard & Silk 2000; Weinberg and Katz 2001). 
If anything, the large fraction of chaotic trajec- 
tories found here and the torques resulting from 
the interaction of two non-axisymmetric struc- 
tures would enhance such coupling — and that 
would include the effect of bar breaking investi- 
gated by Debattista & Sellwood (2000). 

This, however, raises an important question as 
to whether bars would form in the first place, or 
survive long enough, in unmodified CDM halos 
to produce these effects. Athanassoula & Misiri- 
otis (2002) found that not only bars can form 
in centrally-concentrated, but axisymmetric, ha- 
los, but that they are actually much more pro- 
nounced that ones in less concentrated halos. The 
results presented here suggest that small non- 
axisymmetric perturbations may act as to destroy 
these bars. Indeed, we find strong vertical instabil- 
ities at the outer parts of our axisymmetric model 
with centrally-concentrated halo (Model 7) which 
could explain the "X" structure of modeled bars 
when seen edge on. These vertical instabilities 
move inwards and become more pronounced when 
a non-axisymmetric perturbation is present. 

The vertical orbital instability observed was 
found to proceed in tandem with the orbital dis- 
solution in the xy-plane. If the latter is rapid, the 
vertical instability will saturate and is not likely to 
play an important role — since it almost invariably 
takes longer time to develop than the instability 
in the plane. This is especially true if the initial 
distribution is very thin, in which case the vertical 
instability timescale can be very long compared to 
the one in the plane. This is probably the reason 
why Ideta & Hozumi (2000) concluded that the 
vertical instability of an iV-body bar immersed in 
a centrally-concentrated halo does not make an es- 
sential contribution to its evolution. We find that, 
in this case, the trajectories would quickly mix in 
the xy-plane, evolving towards a spatial distribu- 
tion which is isotropic there and forming a lens- 
type configuration. On the other hand, if there 
is sufficient time for the vertical instability of tra- 
jectories to develop, it can play a central role in 
the evolution. For it is especially in the interme- 
diate cases where the instability in the plane may 
not be sufficient in dissolving the bar over short 
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timescales that the additional vertical instability 
may be essential (Section 6.3). The final shape of 
the dissolved object may then be closer in appear- 
ance to a spheroidal component. 

If collective instability, leading to bar forma- 
tion, is active in the case of centrally-concentrated 
triaxial halos, a bar-like structure may appear. 
But if the evolution time of the density distribu- 
tion away from the barred configuration is simply 
related to the Liapunov timescale of the chaotic 
trajectories filling the phase space, such a configu- 
ration will have to be washed out in a few dynam- 
ical times. This should almost certainly be the 
situation in systems where the non-axisymmetric 
non-rotating perturbation reaches 10%. In this 
case there are no regular trajectories supporting 
the bar. Most initial conditions lie in what ap- 
pears to be a connected chaotic region of phase 
space composed of unstable trajectories with e- 
folding time of the order of 10 8 Myr or smaller (cf., 
Fig 14, bottom panel). It is also probably the case 
for systems with a perturbation of 5% (Fig. 14, 
second panel from bottom) . Can the presence of a 
centrally-concentrated triaxial halos in the initial 
stages of galaxy evolution be the cause of appar- 
ent decrease in the bar fraction at redshifts greater 
than 0.5 (e.g., Abraham ct al. 1999)? 

In the case of the "mixed" phase space of sys- 
tems with with smaller non-rotating perturbations 
the situation is more ambiguous. Nevertheless, 
even a system where the non-rotating component 
has a potential axis ratio as large as 0.99 displays 
large connected regions of initial conditions lead- 
ing to strongly chaotic trajectories. A compari- 
son between the bottom panel of Fig. 6 and mid- 
dle panel of Fig. 7 reveals that the configuration 
space structure and orbital stability properties of 
such a model is significantly different from that 
of the model successfully built in a self-consistent 
manner by Pfenniger (1984a, b) — with the latter 
containing far more regular trajectories occupying 
a small configuration space volume aligned with 
the bar. In particular, while both models con- 
tain, at most radii, orbits that are aligned with 
the bar, the variety of available such trajectories 
is far larger in Pfenniger's model. Indeed, simple 
inspection of the spatial structure of a sample of 
bar aligned orbits of Model 6 showed that a large 
fraction of these are too round to contribute to- 
wards a self consistent model. 



The preponderance of chaotic trajectories occu- 
pying large configuration space volumes, and the 
apparent absence of a sufficient population of reg- 
ular bar-supporting orbits, suggest that even mod- 
els with slowly rotating non-axisymmetric pertur- 
bations of ~ 1% cannot have time-independent 
equilibria. The question arises as to whether the 
resulting time dependence would take place on a 
physically relevant timescale. This is not a trivial 
question, since in systems with such mixed phase 
space diffusion times of chaotic trajectories can 
sometimes be very long. However, in the cases 
studied here, unless the isolated regions of regular 
orbits aligned with the bar constitute a sufficient 
contribution towards a self consistent model, the 
bar would have to dissolve in a timescale signifi- 
cantly smaller than a Hubble time. Since we know 
from Section 6.3 that for mildly triaxial centrally- 
concentrated halos the configuration space volume 
of most trajectories is also far larger than the bar's 
over this period. 

It is possible, on the other hand, that central 
cusps in dark halos have been destroyed during 
the formation of the baryonic component via pro- 
cesses involving dynamical friction in inhomoge- 
neous baryon background (El-Zant, Shlosman & 
Hoffman 2001; El-Zant et al. 2002) — or that 
some fundamental physical process relating to the 
nature of the dark matter prevents the formation 
of central density cusps. 
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A. Evaluation of the Liapunov exponents 
A.l. Obtaining the maximal exponent 

Even though Liapunov exponents are a widely known classical tool in nonlinear dynamics, only a few 
studies have used them to quantify chaotic behavior in realistic galactic systems (e.g., Udry & Pfenniger 
1998; Merritt & Fridman 1996; Siopis & Kandrup 2000). There is also the recent study by Fux (2001) on 
the stability of trajectories in a detailed model of the Galaxy. None give a self-contained description of how 
they are obtained. We, therefore, provide such a description here, since these tools are central to the results 
in this paper and in the hope that they would find a wider use in the field. 

Let a system be described by the first order equations of motion (Newtonian second order equations can 
be replaced by two first order ones) 

X = F(X,t) (Al) 

and their variational (linearised) counterparts 

£ = SF. (A2) 

Along a particular trajectory X = X(t, t a , X Q ) against which we would like to measure the deviation, with 
X the initial conditions, (A2) can be rewritten as 

£ = D x F(X(i,i ,X )),£ (A3) 

where D X F is the Jacobian 6N x 6iV matrix dFi/dxj and i,j = 1, 6N. Now let 

X s =X s (X(i,t ,X )) (A4) 

be the fundamental solution of this matrix with the initial condition being the identity matrix. The solution 
of (A3) is then given by (Wiggins 1991) 

€ = X.(i)£ , (A5) 

which describes the evolution under the linearised dynamics with initial conditions £o hi the space of linear 
variations. 

A Liapunov exponent is the infinite limit of the "time-dependent Liapunov exponent" (Wiggins 1991) at 
X in the direction | at time t which is given by 

The Liapunov exponents are then defined as 

<7(£ ,X„)= lim A(€ ,Xo,t). (A7) 

t — >oc 

Numerically of course only A can be calculated. We will refer to the inverse of this time-dependent Liapunov 
exponent as the "exponentiation time," the "exponential timescale" or the e-folding time. 

For a Hamiltonian system with / degrees of freedom, there are 2/ linearly-independent directions in 
phase space for the vector £ to point at, hence there are 2/ Liapunov exponents. A positive Liapunov 
exponent indicates unstable behavior characteristic of chaotic motion. Thus, determining the maximal 
exponent is sufficient for detecting the presence of such behaviour. The evaluation of the maximal exponent 
is straightforward enough. This is because exponential instability, if it is present, will cause almost all 
initial linear tangent space vectors to realign themselves along the subspace of maximal expansion. A 
numerical determination of a Liapunov exponent from almost any initial chosen direction for the linear 
variations will thus tend to give an evaluation of the maximal exponent (Wolf et al. 1985). The only 
complication that arises is that, when the exponentially increasing solutions of the linearised equations 
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become too large, the calculation is slowed down (eventually leading to a numerical overflow) . This is easily 
remedied, however, by application of the "standard algorithm" of Benettin et al. (1976). This algorithm is 
based on the local averaging of the deviation between neighboring states, which is done by dividing the time 
we run the system into n subintervals. An initial linearised deviation £o will, therefore, be transformed into 
X^o, X^Xj£ , X^...X^X^ at times t\, t<2, ■ t n . At iteration n, the part under the logarithm on 
the right hand side of (A6) can then be rewritten as: 



We now successively define 
with 

This means that 



|| X:...X*X 2 s X^ || / || £ || . (A8) 
& = =|| &_i || (A9) 

£i_i = 6-i/ || II- (Aio) 



ii x n ...x 2 x^ n= nr=i ii & ii (aii) 

and, therefore, if one assumes constant intervals At, 

1 i=n log || | 

^standard = -T- hm — . (A12) 

i=i 

In practice, this procedure consists of renormalizing the linearized vector to unity at intervals At, adding 
the logarithm of its norm to the pre-existing sum and restarting the integration with this renormalized unit 
vector serving as initial condition for the variational (linearised) equations. This avoids numerical blowup. 

A. 2. Obtaining the full set 

The problem of the collapse of of the linearized vectors towards the direction of maximum rate of expansion 
can be solved by re-orthogonalizing the vectors at intervals small enough so that linear independence of the 
set of vectors is not completely lost. This can be done by repeated application of the Gramm-Schmidt 
orthogonalization procedure. Therefore, instead of just renormalizing one vector, as in when calculating the 
maximal exponent, at every step of the procedure described above, one re-orthonormalizes a basis set of 
linearly independent tangent space vectors (£*,£?, to obtain a new set of orthonormal vectors given 

by 

ii = Sty, ( A13 ) 



i\ 



,2 [ii 'ii I ii 

ii = . rrr : ( A14 ) 



e-(^r)g--.-(^.i;)i- 
• " iu;'-(«s'.r i )«r , -...-(^.«, , )«'»' 

That is the new set of vectors is made orthonormal by simply normalizing the first vector, then subtracting 
the projection of the first vector on the second and normalizing to get the new second vector. After that, 
the first two new vectors are subtracted form the third vector of the original set which is then normalized to 
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obtain the new third vector, etc. The vector ^ i continues to seek out the direction of maximum expansion 

(since it has the same direction as £ while ^ ) , and $ ti span the most rapidly growing two-dimensional 
subspace, and, in general, the first k orthogonalized vectors span the same subspace as the first k vectors of 
the original set. Also since the new set of vectors is orthogonal, one may determine the Liapunov exponents 
from the mean rate of growth of the projection of the new vectors on the old ones. A FORTRAN routine 
that finds the Liapunov exponents using this procedure is given by Wolf et al. (1985). 

Better numerical stability can be achieved by using a modified Gramm-Schmidt orthogonalization proce- 
dure (e.g., Lawson & Hanson 1974; Noble 1988). This amounts to the following: instead of subtracting the 
k — 1 preceding vectors from the k th vector, thus making the latter orthogonal to all of the former, one starts 
with the k th vector and makes all the following 2f — k vectors orthogonal to that vector. Thus we start with 
the first vector to get 

il 
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& = (A18) 
This is done for every vector until we get to 

& = ^V" ( A1 9) 
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Here we use only this more stable modified algorithm. 



B. The Kolmogorov entropy and statistical evolution 

A discussion of the precise conditions under which the Liapunov exponents are good indicators of statis- 
tical behavior can be found elsewhere (e.g., Pesin 1989; Eckmann & Ruelle 1985). We just mention that 
under certain conditions believed to be satisfied for many physical systems the exponents are related to the 
Kolmogorov-Sinai entropy by 

KS= f X )dX (Bl) 

•* i 

where the sum is taken over all positive exponents and the integral is over all possible initial conditions. For 
a single trajectory the KS entropy (in this case simply the sum of the positive exponents) is a measure of the 
information loss about the initial phase space point X as the trajectory propagates, or the "complexity" 
of the trajectory. It is zero for regular orbits, where the motion is separated into recurring oscillations in 
each degree of freedom. For a set of trajectories in a connected chaotic region of phase space, it is a measure 
of the rate of increase of coarse-grained phase volume they represent, and hence increase in the statistical 
entropy. It is, therefore, possible to interpret Liapunov exponents as measuring the rate of evolution of an 
initially improbable distribution of phase space points. Nevertheless because of the existence of phase space 
barriers in systems with mixed phase spaces — that is ones with coexisting regular and chaotic regions — 
some trajectories may take very long to reach an invariant distribution. The correspondence is therefore not 
always so straightforward. It was one of our goals in this paper to test it for our current systems. It turns 
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out that, for these systems, the exponents are very useful diagnostics of (at least) the configuration space 
structure of trajectories over timescale of the order of a Hubble time.. 
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